How multiple air pollutants affect hand, foot, and mouth disease incidence in children: assessing effect modification by geographical context in multicity of Sichuan, southwest China

Background Several studies have suggested a significant association of hand, foot, and mouth disease (HFMD) with ambient air pollutants. Existing studies have characterized the role of air pollutants on HFMD using only risk ratio measures while ignoring the attributable burden. And whether the geographical context (i.e., diverse topographic features) could modulate the relationships is unclear. Methods Daily reported childhood HFMD counts, ambient air pollution, and meteorological data during 2015–2017 were collected for each of 21 cities in Sichuan Province. A multistage analysis was carried out in different populations based on geographical context to assess effect modification by topographic conditions. We first constructed a distributed lag nonlinear model (DLNM) for each city to describe the relationships with risk ratio measures. Then, we applied a multivariate meta-regression to estimate the pooled effects of multiple air pollutants on HFMD from the exposure and lagged dimensions. Finally, attributable risks measures were calculated to quantify HFMD burden by air pollution. Results Based on 207554 HFMD cases in Sichuan Province, significant associations of HFMD with ambient air pollutants were observed mainly at relatively high exposure ranges. The effects of ambient air pollutants on HFMD are most pronounced on lag0 or around lag7, with relative risks gradually approaching the reference line thereafter. The attributable risks of O3 were much greater than those of other air pollutants, particularly in basin and mountain regions. Conclusions This study revealed significant pooled relationships between multiple air pollutants and HFMD incidence from both exposure and lag dimensions. However, the specific effects, including RRs and ARs, differ depending on the air pollution variable and geographical context. These findings provide local authorities with more evidence to determine key air pollutants and regions for devising and implementing targeted interventions. Supplementary Information The online version contains supplementary material available at 10.1186/s12889-023-17484-9.


Background
Hand, foot, and mouth disease, an acute intestinal infectious disease, is prevalent among children below 5 years old worldwide [1].During recent decades, multiple HFMD epidemics have occurred mostly in the Southeast regions, with considerable numbers of infections.Most HFMD cases are self-limiting and characterized by common symptoms such as fever and rash, while a minority of infected children, especially those under the age of 3, may experience severe complications of the central nervous system and circulatory system, which can even lead to death [2].A variety of human enteroviruses can cause HFMD through respiratory droplets, the fecal-oral route, and close contact.Human enterovirus A71 (EV-71) and coxsackievirus A16 (CVA16) are the primary pathogenic serotypes that lead to HFMD outbreaks [3].However, in recent years, the dominant causative agent of HFMD has undergone a shift with the rapid surge in cases caused by CV-A6 and CV-A10 [4,5].Regrettably, effective antiviral drugs for the prevention and treatment of HFMD are still unavailable.In 2016, several monovalent EV-A71 vaccines were launched in China.These vaccines exhibited high effectiveness in preventing HFMD associated with EV-A71 and mitigating the severity of HFMD in infected patients [6,7].Nonetheless, individuals susceptible to HFMD must pay the full price of the vaccines to acquire protection exclusively against EV-A71 without efficient cross-protection against other enterovirus serotypes [6,8].Even after the approval of EV-A71 vaccines, the HFMD remains a high prevalent disease which poses an important public health threat to child health [9][10][11].HFMD persists with a much higher incidence rate than other notifiable infectious diseases in mainland China, with more than one million infected cases each year [9][10][11][12][13].
HFMD exhibits distinct seasonal patterns [14].Previous epidemiological studies have demonstrated the relationship between environmental variables, such as temperature, relative humidity, and other meteorological factors, and HFMD incidence [15,16].In recent years, with the rapid development of economic growth, urbanization, and industrialization, the effects of ambient air pollutants on health have gained considerable attention from researchers [17,18].Exploring the impacts of air pollutants on HFMD can help clarify the potential factors, which further provide environmental evidence for targeted policy recommendations and prevention strategies.Nonetheless, several studies have quantified the associations of HFMD with multiple air pollutants, which showed regional inconsistency.For example, a study in Hefei revealed a linear relationship between SO 2 and HFMD, with the highest cumulative effect at lag day 0-1 among scattered children [19].In contrast, Gu et al.proposed a heterogeneous J-shape relationship between HFMD and SO 2 in Ningbo [20].In addition, O 3 has been found to have an inverse V-shape relationship with HFMD in Ningbo [20], while an M-shape pattern was observed in Shenzhen [21].In addition, studies were usually carried out in a single city with moderate or mild air pollution and applying methods with different relationship assumptions or model specifications.Hence, the generalizability and comparability of the findings were limited across multiple locations.Furthermore, previous studies relied solely on risk ratio measures such as RR and OR to evaluate the impact of ambient air pollution variables on HFMD, which provided only partial information about the excess burden of HFMD caused by air pollution.To date, there is no evidence on the excess burden of HFMD attributable to air pollution.Therefore, it is necessary to rely on multiple insights viewing the role of air pollutants in HFMD through appropriate multicity strategies and attributable indices (i.e., AF and AN) to help elucidate the possible influencing factors and supplement the evidence for health effect estimation.
Epidemiological studies have shown that a few cityspecific characteristics, such as socioeconomic status and climatic conditions, could modify the associations of HFMD with ambient air pollution [22,23].The topographic condition of a city is a critical characteristic that can reflect its unique natural environment and have direct and indirect impacts on the social and economic conditions of the city.However, whether topography acts as a comprehensive factor affecting the association of ambient air pollutants with HFMD has not yet been investigated.Sichuan Province, located in Southwest China, encompasses a variety of topographic conditions, including basins, plateaus, and mountains.Meanwhile, Sichuan Province has imbalanced socioeconomic development and diverse climate conditions across different topographies [24].Consequently, the distributions of air pollution concentrations vary greatly across its topographies in Sichuan Province [25].Sichuan Basin, in particular, has emerged as a most heavily air polluted region in China [26].Therefore, Sichuan Province could serve as a representative region to investigate and compare the complex associations between air pollution and HFMD across different geographical areas under the combined impacts of the environment and social economy, particularly to provide novel quantitative evidence in heavily polluted areas.
This study aimed to explore and compare the relationships in the context of heavy air pollution utilizing DLNM in three geographical regions (i.e., Sichuan Basin, West Sichuan Plateau, and Southwest Sichuan Mountain): a) the pooled cumulative effects of air pollution on HFMD; b) the pooled lag effects of air pollution on HFMD; c) the attributable risk of air pollution on HFMD.

Setting and data
Sichuan Province, composed of 21 cities, covers a total area of more than 486,000 square kilometers in Southwest China.Sichuan Province has a diverse and complex topography with five typical landforms, including basins, plains, hills, mountains, and plateaus.Additionally, Sichuan Province exhibits persistent regional variations in climate.The eastern Sichuan Basin is located in the subtropical humid climate zone.The southwestern Mountain and western Plateau are located in the subtropical and high-altitude cold-climate zones, respectively [27].Considering the unique climatic and geographical characteristics, the national guide of divisions published by the China Meteorological Administration [28] divides Sichuan Province into 3 regions: the Sichuan Basin (Chengdu, Zigong, Deyang, and 14 other cities), West Sichuan Plateau (Aba and Ganzi), and Southwest Sichuan Mountain (Liangshan and Panzhihua) (Fig. 1).

Fig. 1 Geographic locations for the 21 cities and air quality monitoring stations in Sichuan Province
On May 2, 2008, HFMD was cataloged into a notifiable infectious disease classified C under the unified report and standard management.The surveillance data for daily HFMD cases during 2015-2017 were retrieved from the national information system for disease control and prevention managed by the China CDC.Each reported patient can be extracted following information, including demographic value (gender, date of birth, and address), case classification, and date of onset.We compiled all of the case information into daily time series data of HFMD counts for each of the 21 cities in Sichuan Province.As previous studies demonstrated that children below 5 years old were the most vulnerable group to infect HFMD, we only included patients aged 0-5, which accounted for an overwhelming majority of HFMD cases (> 96%) according to the preliminary analysis.
Daily air pollution data, including PM 10 (μg/m 3 ), PM 2.5 (μg/m 3 ), SO 2 (μg/m 3 ), NO 2 (μg/m 3 ), O 3 (μg/m 3 ), CO (mg/ m 3 ), and air quality index (AQI), for the same period in 21 cities were downloaded from the Sichuan Environmental Monitoring Center.The monitoring data are of high quality with a proportion of missing data less than 1%.Linear interpolation was used to fill in missing values.
Daily meteorological data from 41 monitoring stations in Sichuan Province were provided by the official website of National Meteorological Data Sharing Center.Seven meteorological indicators, namely, mean temperature (°C), minimum and maximum temperature (°C), relative humidity (%), air pressure (hPa), sunshine duration (h), and wind velocity (m/s), were considered in this study.The proportion of valid data for each station is greater than 99%.To fill the missing values, we used different methods according to the nature of each meteorological indicator.Missing values of sunshine duration were replaced with zero.For the other six variables, we used linear interpolation to replace missing values.To estimate daily meteorological exposures at the city level in Sichuan Province, we implemented the kriging interpolation method, allowing us to make full use of the information from all meteorological stations.
Ultimately, daily series of HFMD counts, air pollution concentrations, and meteorological data of each of the 21 cities in Sichuan Province were matched through unified city-specific codes for further analysis.

Statistical analysis
In this study, a multistage and multicity time series regression was conducted to describe the effects of ambient air pollution variables on HFMD from multiple perspectives.In the first stage, the estimated effects of air pollution from two dimensions (i.e.exposure-response and lag-response) were quantified for each of the 21 cities in Sichuan Province.In the second stage, pooled exposure-response and lag-response relationships at the regional level were described by a multivariate metaanalysis based on the city-specific estimates from the first stage [29,30].In the last stage, the attributable risks of HFMD associated with air pollution were calculated based on the backward perspective [31].

DLNM construction
Due to the frequent occurrence of sporadic HFMD cases or outbreaks, the count data were overdispersion.Therefore, we employed the distributed lag nonlinear model assumed a quasi-Poisson distribution to quantify the impacts of multiple air pollutants on HFMD.According to our initial analysis, we found a strong correlation between PM 10 , PM 2.5 , and AQI, which could affect the accuracy of the results.To improve the comparability with previous studies and consider the possible hypotheses on biological mechanisms, only PM 10 was retained in this study for further analysis among these highly correlated variables.We established uniform model structures and parameters for each city by considering prior knowledge and conducting a series of sensitivity analyses based on the multi-pollutant model (see more details in Additional file).The final model for each city is expressed as follows: where α is the intercept; i is the city code (i = 1,2,3,…,21); Y it represents the number of daily HFMD counts in city i corresponding to day t (t = 1,2,3,…,1096); Quasi-Poisson refers to the quasi-Poisson model constructed; and cb(air pollutants it ) is the cross basis function of PM 10 , SO 2 , NO 2 , CO, and O 3 associated with HFMD.Specifically, the two dimensions (i.e.exposure-response and lag-response) were described by the natural cubic spline (ns) with the degrees of freedom (df ) setting to 3, and the lag period was set to 0-14 days to capture the lag effects of air pollution.We defined a natural cubic spline using 8 dfs per year expressed as ns(day,8*3) in the model to remove the confounding of seasonality and long-term trend.A nonlinear term of temperature was incorporated by first summarizing multiple lags (0-14) using a simple moving average (SMA) and then applying a ns with 3 dfs

Multivariate meta-regression
The exposure-lag-response associations of HFMD with air pollution variables were estimated in the first stage and then reduced to unidimensional vectors for the following multivariate meta-analysis.For each air pollutant, the lag-response relationships at specific air pollution levels, including the 25 th and 75 th percentiles, were extracted.In addition, the delayed effects at each air pollution level were summed to represent the cumulative air pollutants-HFMD relationship.Taking the region as the unit, a intercept model of the multivariate metaregression was used to meta-analyze the city-specific effect estimates after dimensionality reduction to obtain the region-specific pooled lag-response relationship and cumulative exposure-response relationship.

Attributable risks calculation
The air pollution with the lowest RR, derived from the city-specific pooled air pollutant-HFMD curve ranging from the 5 th to 95 th percentile of the air pollutant, was taken as the counterfactual condition to calculate the attributable risk [32,33].For each day of the series, the cumulative RR corresponding to the air pollutant concentration obtained from the first stage was used to calculate the attributable HFMD cases and fraction in each city [31].To quantify and compare the impacts of each air pollutant at different exposure levels, the 5 th , 10 th , 25 th , 50 th , 75 th , 90 th , and 95 th percentiles were selected as the cutoff points to calculate the ARs in low exposure intervals (P 0 ~P5, P 0 ~P10, P 0 ~P25, P 0 ~P50 ) and high exposure intervals (P 50 ~P100, P 75 ~P100, P 90 ~P100, P 95 ~P100 ) [34].
The numbers of HFMD cases attributable to high and low exposure levels of air pollution were calculated by summing the components for all of the days in the time series that met the requirements.The region-specific ANs were further calculated and compared by summing the subsets corresponding to cities located in Sichuan Basin, West Sichuan Plateau, or Southwest Sichuan Mountain.The attributable fraction was produced by dividing the AN by the corresponding total number of HFMD cases.Monte Carlo simulations were utilized to obtain the empirical confidence intervals for AFs and ANs under the assumption of multivariate normal distribution.All statistical analyses were performed in R 3.5.1.We adopted R packages "dlnm" and "mvmeta" to fit the DLNMs and multivariate meta-regression models, respectively.Two-tailed P < 0.05 were defined as statistically significant for all statistical tests.

Descriptive statistics of research data
The dataset included 207554 cases under 5 years old with HFMD in Sichuan Province during 2015-2017, with yearly totals of 62613, 86864, and 58077 cases, respectively.The HFMD cases were mostly clustered in the eastern basin, which has a high population density and relatively developed economy.We found a similar seasonal distribution of HFMD in both the Sichuan Basin and the West Sichuan Plateau.Two peaks were observed during spring (from April to June) and early winter (from November to December) (Fig. 2).In 2016, both the annual total number and the peak number of cases were higher than those in the other two years.Furthermore, we found that PM 2.5 , PM 10 , and AQI had similar distribution patterns across the different regions.High concentrations of particulate matter typically accumulated in the Sichuan Basin, while the average concentrations of SO 2 , CO, NO 2 , and O 3 were higher in Southwest Sichuan Mountain (Table 1).

Pooled effects of ambient air pollution variables on HFMD in the exposure dimension
Our results showed significant associations of air pollutants on HFMD in three geographical divisions.In Sichuan Basin, we found cutoff values of 38.4 μg/m 3 and 58.3 μg/m 3 for the association of HFMD with SO 2 and NO 2 , which correspond to the slopes of the curves nearest to zero.After that, the cumulative relative risk increased approximately linearly with SO 2 and NO 2 .In contrast, PM 10 displayed an approximately linear negative correlation with pediatric HFMD when the concentration exceeded the median.Similarly, in West Sichuan Plateau, once CO exceeded the median, the relative risk of HFMD began to decrease linearly.In Southwest Sichuan Mountain, we observed a significant nonlinear association between O 3 and HFMD.The risk decreased with increasing O 3 initially, reaching the lowest RR of 0.412 (95% CI: 0.224-0.759)at 127.8 μg/m 3 of O 3 , and then the curve began to flatten (Fig. 3).In summary, significant associations were observed among all the air pollutants with HFMD in three geographical divisions, but the detailed shape patterns for the air pollutant-HFMD relationships differed across geographical divisions.

Comparison of the attributable risks of air pollutants on HFMD
Initially, we compared the attributable fractions for air pollutants at different concentration ranges in three geographical divisions (Fig. 5).We found that significant AFs of PM 10 , SO 2 , and NO 2 mostly in Sichuan Basin.The health effects of O 3 were much greater than those of other air pollutants, particularly in basin and mountain regions.Detailedly, the estimated AFs at the low-level PM 10 concentration intervals ranged from 1.04% (95% CI: 0.24%-1.69%)to 15.88% (95% CI: 9.3%-20.72%) in the Sichuan Basin.Likewise, the AFs for the low-level CO concentration intervals were 2.38% (95% CI: 1.19%-3.43%),5.10% (95% CI: 2.65%-7.14%),9.19% (95% CI: 5.17%-12.46%),and 15.55% (95% CI: 8.48%-20.21%),higher than those for the   5).Due to the large population and a huge number of HFMD infections in Sichuan Basin, the ANs in this region for each air pollutant were also much higher than those in West Sichuan Plateau and Southwest Sichuan Mountain (Table 2 and 3).In Sichuan Basin, the ANs for high-level SO 2 and NO 2 account for most of the HFMD burden compared with other air pollutants.Compared with the counterfactual condition, the extreme highlevel PM 10 concentration range showed a protected effect from HFMD with an AN of 1778 (Table 2).In West Sichuan Plateau, significant ANs were only observed for high-level NO 2 and CO concentration intervals above the 50 th and 75 th percentiles for low-level NO 2 and CO concentration intervals below the 10 th percentile.In Southwest Sichuan Mountain, O 3 was the most remarkable air pollutant with the highest ANs for low-level O 3 concentration intervals ranging from 550 (95% CI: 307-741) to 3011 (95% CI: 1786-3747).

Discussion
In recent years, air pollution, one of the key environmental problems affecting human health, has received wide attention from researchers.Air pollutants can increase the risk of morbidity of infectious and noninfectious diseases in the respiratory, circulatory, and nervous systems [35][36][37].In this study, a multistage strategy based on DLNMs was applied to comprehensively explore the possible nonlinear and delayed associations and quantify the attributable risks between air pollutants and HFMD.The results suggested significant pooled nonlinear associations of PM 10 , CO, and O 3 with HFMD and approximately linear relationships of SO 2 and NO 2 with HFMD.The pooled effects of air pollution variables on HFMD are most pronounced on lag0 or around lag7.These Fig. 3 Significant overall pooled relationships pattern between multiple air pollutants and HFMD at the regional level effects, including RRs and ARs, were closely related to geographical divisions and exposure levels.To the best of our knowledge, this study is the first to characterize the pooled lag effects of air pollution at the regional level and quantify the burden of air pollution on HFMD by calculating attributable risks.
This study suggests that further attention should be given to the effects of air pollutants in areas with high exposure concentrations.It was observed that the overall effects of PM 10 , SO 2 , and NO 2 on HFMD were significant only in the basin area, while those of CO and O 3 on HFMD were significant only in the plateau and mountain areas, respectively.Specifically, the risk of HFMD gets higher with increasing concentrations of SO 2 and NO 2 , which shows consistency with previous studies.For instance, a study in Ningbo suggested that high concentrations of SO 2 and NO 2 elevated the risk of HFMD [20].Shao et al. indicated that HFMD incidence is related to cytokines (e.g., TNF-α, IFN-β, IL-4, IL-12, IL-18) and chemokines, which are usually affected by SO 2 [38].As a result, high SO 2 remarkably weakens the body resistance of children to viral infection.Particularly under the typical topography of the Sichuan Basin with high humidity throughout the year, due to the high solubility of SO 2 in water, inhaled SO 2 and its in vivo derivatives may rapidly enter the bloodstream and produce toxic effects in the respiratory system, decreasing the host defense responses to enteroviruses [39,40].Many studies have shown that NO 2 influences the incidence of respiratory diseases in children, such as low lung function in school children [41] and childhood asthma [42,43].However, research on the biological mechanism and epidemiological relationship between NO 2 and HFMD remains limited.As one of the key components of air pollution, NO 2 is closely related to vehicle exhaust emissions.With the rapid pace of urbanization, car ownership in all cities has Fig. 4 The regional pooled lagged effects of air pollutants on HFMD in three geographical regions expanded dramatically, especially in the Sichuan Basin, which leads in both level and speed of development.The average concentration of NO 2 in the basin is 30.5 μg/ m 3 , which is much higher than that in other regions in Sichuan Province.Therefore, in the basin area, it is necessary to comprehensively consider the exposure concentrations of NO 2 and SO 2 , along with the specific curve pattern of HFMD with air pollutants, to formulate targeted health guidance strategies.
This study found that high concentrations of PM 10 , O 3 , and CO were related to a lower risk of HFMD, with each significant air pollutant being identified in a specific geographical region.It is believed that particulate matter may facilitate the adhesion and transmission of enterovirus and that it may also influence the inflammatory response, leading to increased susceptibility to infection [44,45].However, due to poor visibility on polluted days with high concentrations of particulate matter, outdoor activities for children may be limited to reduce the risk of contact with enterovirus and particulate matter in the environment, ultimately leading to a reduction in the risk of HFMD associated with high particulate matter [46,47].Some studies conducted in Ningbo [48] and Shenzhen [21] found no significant relationship between PM 10 and HFMD, which differs from the results of this study.This may be due to the fact that Sichuan Basin is a typical region with a high incidence of HFMD and much higher concentrations of particulate matter than the coastal cities, making the effect of PM 10 more likely to be identified.Furthermore, evidence from another heavily polluted city, Shijiazhuang, suggests that high particulate matter may have a protective effect [49].Therefore, we can reasonably speculate that the average concentrations of air pollutants in cities might be a potential factor that could impact the relationship pattern of HFMD with air pollution.Similarly, a significant negative association of O 3 with HFMD was found in the mountainous area with high O 3 concentrations.Ozone can inhibit the survival and replication of bacteria and viruses in Fig. 5 Attributable fractions of HFMD incidence by multiple air pollutants at different concentration ranges in three geographical divisions in Sichuan Province the external environment.Some studies found that the inactivation of enterovirus EV-A71 affected by O 3 is related to the dissolution dynamics of O 3 [50,51].In contrast, the protective effect of high CO was identified in the plateau area with lower concentrations than other regions.Some studies have indicated that chronic exposure to CO may increase the risk of cardiovascular disease [52], while others have suggested that shortterm exposure to low-level CO presents a protective effect and reduces the risk of hospital admission for  respiratory infectious diseases [53].The discrepancies among these results might be related to the differences in disease outcomes, geographic locations, and socioeconomic conditions.Further evidence is still needed to detailedly explore the potential reasons and better characterize the relationship between CO and HFMD.Our findings reveal that the effects of ambient air pollution indicators on HFMD are most pronounced on lag0 or around lag7, with relative risks gradually approaching the reference line thereafter.Utilizing the lag distribution of the effects could prove advantageous in directing the prevention and control of HFMD.It is noteworthy that the effect of O 3 was observed on the third lag day and lasted for over a week.However, single-day lag effects provided limited information.To quantify and compare the effects of air pollutants, we further calculated attributable fractions and attributable numbers under various exposure ranges and observed that the attributable risks of PM 10 , SO 2 , and NO 2 tended to be higher in the Sichuan Basin.This may be related to the differences in the distribution of air pollutants.A comparison of the attributable fractions of air pollutants revealed that the health effects of O 3 were much greater than those of other air pollutants, particularly in basin and mountain regions.In previous studies, temperature and relative humidity, the two most popular environmental factors, were commonly identified to have significant impacts on HFMD risks [33,54,55].The attributable fractions of high temperature and relative humidity in previous studies reached 39.55% [54] and 14.56% [33], respectively.The latter was relatively close to the AFs for the high levels of PM 10 (13.38%) and O 3 (16.83%)that were calculated in this study.More attention should be given to the associations between HFMD and air pollutants, especially PM 10 and O 3 , next only to temperature.In summary, in addition to the bidimensional pooled exposure-lag-response relationships, the measurement of attributable risk can provide new insights of association from a quantitative perspective, particularly with advantages in the comparison of HFMD burden attributed to different air pollutants.
Several limitations should be noted in this study.First, the intrinsic nature of this ecological study precludes the establishment of a causal relationship between ambient air pollution and HFMD, despite the identification of short-term associations and attributable risks.Nonetheless, the results may serve as the basis for future research hypotheses.Second, to leverage the air pollution variables to their fullest extent, data were collected only after 2015.Despite this, this study still possessed sufficient test power given the daily data that were utilized.Third, some measurement errors or inaccuracies such as under-reporting HFMD incidence may have been present in the passive surveillance data, which potentially leads to underestimation of the effects [56].Fourth, only partial available and measurable meteorological and air pollution variables were considered as confounders in the models.There may have been additional potential confounders or external factors that are worth considering in further research.

Conclusions
This study revealed a significant two-dimension relationship (i.e.exposure-response and lag-response) between several air pollution indicators and HFMD incidence.However, the specific effects, including RRs and ARs, differ depending on the air pollution variable and topographic division.Children living in the Sichuan Basin are more vulnerable to PM 10 -related HFMD compared with SO 2 and NO 2 .The health effects of O 3 were much greater than those of other air pollutants, particularly in basin and mountain regions.These findings could furnish local authorities with more evidence to determine key air pollutants and regions for devising and implementing targeted interventions.

Fig. 2
Fig. 2 Daily distribution of HFMD cases from three regions in Sichuan Province, 2015-2017 to account for the confounding of temperature.Other meteorological factors (relative humidity, wind velocity, and sunshine duration) were included in the model by a simple-moving-average variable as SMA(Other confounder).The respective coefficient of the other confounder was denoted as δ oc .We used the logarithm scale of HFMD counts at lag1 and lag2 as the autoregressive terms to adjust the residual autocorrelation.The day of week (dow t ) and holidays (holiday t ) were included as dummy variables.And η, β, and γ are coefficients of corresponding terms.Childpop im was an offset term for the child population of the city i on year m.The median values of PM 10 , SO 2 , NO 2 , CO, and O 3 across all the cities (i.e., 60.0 μg/m 3 , 13.0 μg/m 3 , 27.0 μg/m 3 , 0.8 mg/m 3 , 57.0 μg/m 3 ) were defined as the reference to calculate relative risks.

Table 1
Descriptive statistics of air pollution variables and HFMD cases from three regions in Sichuan Province, 2015-2017

Table 2
Numbers of HFMD incidence attributable to multiple air pollutants at high concentration ranges in three geographical divisions in Sichuan Province

Table 3
Numbers of HFMD incidence attributable to multiple air pollutants at low concentration ranges in three geographical divisions in Sichuan Province